library(tidyverse)
library(cowplot)
library(lubridate)
library(mgcv)
source("UVP_2017_library.R")
theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
bes<- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
bds <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
ues <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m,
depth = [32mcol_double()[39m,
psd_gam = [32mcol_double()[39m,
vol = [32mcol_double()[39m,
sizeclass = [31mcol_character()[39m,
lb = [32mcol_double()[39m,
ub = [32mcol_double()[39m,
binsize = [32mcol_double()[39m,
TotalParticles = [32mcol_double()[39m,
nparticles = [32mcol_double()[39m,
n_nparticles = [32mcol_double()[39m,
biovolume = [32mcol_double()[39m,
speed = [32mcol_double()[39m,
flux = [32mcol_double()[39m,
flux_fit = [32mcol_double()[39m,
GamPredictTP = [32mcol_double()[39m
)
uds <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
PhoticBase <- 160
OMZBase <- 850
PlotNParticles <- uds %>%
ggplot(aes(x = tot_nparticles, y = depth, col = profile)) +
facet_wrap(~project) +
geom_point(alpha = 0.3, shape = 1) +
scale_y_reverse() + scale_x_log10()
PlotNParticles
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
#FSG4 <- gam(tot_nparticles~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.0025 0.1091 82.55 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.1031 1.196 3.873 0.0583 .
s(Day) 1.5509 1.796 3.715 0.0801 .
s(Hour) 0.4267 2.000 0.281 0.2713
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.134 Deviance explained = 17.9%
GCV = 0.76575 Scale est. = 0.71367 n = 60
#summary(FSG2)
#summary(FSG3)
#summary(FSG4)
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.006778976
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.07324885
summary(FSG3)$r.sq
[1] 0.0542546
But there is between projects:
ProjGam <- gam(tot_nparticles~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.9983 0.4445 20.24 <2e-16 ***
factor(project)P16 17.0001 1.3754 12.36 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1 1 3.552 0.064 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.7 Deviance explained = 70.9%
GCV = 12.411 Scale est. = 11.855 n = 67
library(scales)
#https://stackoverflow.com/questions/30179442/plotting-minor-breaks-on-a-log-scale-with-ggplot
log_breaks = function(maj, radix=10) {
function(x) {
minx = floor(min(logb(x,radix), na.rm=T)) - 1
maxx = ceiling(max(logb(x,radix), na.rm=T)) + 1
n_major = maxx - minx + 1
major_breaks = seq(minx, maxx, by=1)
if (maj) {
breaks = major_breaks
} else {
steps = logb(1:(radix-1),radix)
breaks = rep(steps, times=n_major) +
rep(major_breaks, each=radix-1)
}
radix^breaks
}
}
scale_x_log_eng = function(..., radix=10) {
scale_x_continuous(...,
trans=log_trans(radix),
breaks=log_breaks(TRUE, radix),
minor_breaks=log_breaks(FALSE, radix))
}
#theme_set(theme_bw)
PlotPSDmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(y = "Depth (m)", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
#theme_set(theme_cowplot)
PlotParticlesmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
scale_x_log10(breaks = c(10, 100, 1000), minor = c(5, 50, 500)) +
#theme(legend.position = "none") +
#scale_x_log_eng()+
labs(y = "Depth (m)", x = "Particles / L") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
# PlotFluxmany <- uds %>%
# filter(project == "ETNP") %>%
# ggplot(aes(x = tot_flux_fit, y = depth, shape = factor(day(time)), fill = hour(time))) +
#
# #geom_path(aes(x = psd_gam)) +
# #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
# geom_point(alpha = .6, size = 2, stroke = 1) +
# scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
# scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
# scale_x_log10() + theme(legend.position = "none")
plot_grid(
PlotParticlesmany,
PlotPSDmany,
rel_widths = c(2, 3)
)
Removed 266 rows containing missing values (geom_point).Removed 266 rows containing missing values (geom_point).
ggsave("figures/ParticlesPSDMany.png")
Saving 12 x 7.41 in image
ggsave("figures/ParticlesPSDMany.svg")
Saving 12 x 7.41 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG4 <- gam(psd~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01975 -200.6 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.719 1.921 56.605 3.83e-15 ***
s(Day) 1.000 1.000 1.278 0.2631
s(Hour) 1.547 2.000 2.763 0.0347 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.639 Deviance explained = 66.5%
GCV = 0.025652 Scale est. = 0.023401 n = 60
#summary(FSG2)
#summary(FSG3)
summary(FSG4)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01978 -200.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.721 1.922 56.019 4.39e-15 ***
s(Hour) 1.609 2.000 3.023 0.0294 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.638 Deviance explained = 65.9%
GCV = 0.025297 Scale est. = 0.023471 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.03439091
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.003377312
summary(FSG3)$r.sq
[1] 0.6015425
Not a significant difference in PSD with respect to time.
But there is between projects:
ProjGam <- gam(psd~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96167 0.02491 -159.029 <2e-16 ***
factor(project)P16 -0.19977 0.07708 -2.592 0.0118 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.224 1.398 34.54 2.26e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.437 Deviance explained = 45.6%
GCV = 0.039117 Scale est. = 0.037234 n = 67
I wonder if I can show that the profiles aren’t statistically significanlty different. Or that they are for that matter… I think in that case, I run a gam with and without a parameter for profile… And then quantify the effect size of that parameter
Or follow this Gavin Simpson Post https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/
or anova.gam {mgcv}
Calculate gams for each profile, and then run anova.gam to see if they are different…
PlotNParticlesEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
#geom_path(aes(x = tot_nparticles)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
labs(x = "Particles/L", y = "Depth (m)") +
theme(legend.position = "none") +
scale_shape_manual(values = c(1:5)) +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
PlotNParticlesEP
I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)
PlotPSDEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
geom_path(aes(x = psd_gam)) +
geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_color_manual(values = c("gray20", "brown")) +
scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
PlotPSDEP
I may just cow these togther.
plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))
Removed 1211 rows containing missing values (geom_point).Removed 1211 rows containing missing values (geom_point).Removed 1211 row(s) containing missing values (geom_path).
ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")
Saving 10 x 4 in image
ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")
Saving 10 x 4 in image
mainParticleComponents <- bds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
select(project, profile, depth,
tot_nparticles, small_nparticles, big_nparticles,
tot_psd = psd, small_psd, big_psd,
tot_flux_fit, small_flux_fit, big_flux_fit) %>%
pivot_longer(cols = -c("project", "profile", "depth")) %>%
separate(name, c("size", "meas")) %>%
mutate(meas = recode(meas, nparticles = "particles/L")) %>%
mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))
Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].
PlotFlx <- mainParticleComponents %>%
filter(meas != "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
PlotPSD <- mainParticleComponents %>%
filter(meas == "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))
Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).
ggsave("figures/BigVsSmall.svg")
Saving 7.29 x 4.5 in image
ggsave("figures/BigVsSmall.png")
Saving 7.29 x 4.5 in image
Flux small and flux tot track so closely because ag > psd. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top. time 2017-01-13 17:51:31 202.4 L
eg_dataline <- bds %>%
filter(profile == "stn_043", depth == 162.5)
eg_slope = eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)
eg_datablock <- bes %>%
filter(profile == "stn_043", depth == 162.5)
eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)
eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)
EgNNP <- eg_datablock %>%
ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")
EgNP <- eg_datablock %>%
ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")
EgTP <- eg_datablock %>%
ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")
plot_grid(EgNNP, EgTP, labels = c("A", "B"))
Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis
ggsave("figures/ExamplePSD163m.png")
Saving 7.29 x 4.5 in image
ggsave("figures/ExamplePSD163m.svg")
Saving 7.29 x 4.5 in image
bds %>%
ggplot(aes(y = depth, x = Flux_Smooth, col = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10()
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(Flux_Smooth~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(Flux_Smooth ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9312 37.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.874 1.984 12.034 0.000111 ***
s(Day) 1.828 1.969 2.774 0.089781 .
s(Hour) 1.515 2.000 2.784 0.031360 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.333 Deviance explained = 39.2%
GCV = 58.038 Scale est. = 52.024 n = 60
summary(FSG2)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9774 35.53 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.864 1.981 10.835 0.000255 ***
s(Day) 1.737 1.931 1.586 0.243284
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.265 Deviance explained = 31%
GCV = 62.081 Scale est. = 57.32 n = 60
summary(FSG3)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9944 34.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.857 1.979 10.73 0.000277 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.239 Deviance explained = 26.3%
GCV = 62.291 Scale est. = 59.325 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.06791706
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.02571191
summary(FSG3)$r.sq
[1] 0.2392424
bds %>% filter(project == "ETNP") %>% select(profile, depth, Flux_Smooth) %>% pivot_wider(names_from = profile, values_from = Flux_Smooth)
Something is off. All of the flux profiles are identical. Skip this
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
plt1 <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = DFP, col = factor(time), shape = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + xlim(c(0.5, 1.5))+ geom_vline(xintercept = 1) +
scale_color_manual(values = c(rep("black", 5), rep("blue", 5))) + scale_shape_manual(values = rep(1:5, 2))
plotly::ggplotly(plt1)
What the heck is going on with DFP here. Why is it usually > 1 shouldn’t it be less than 1 when flux is decreasing? This very deep increasing flux seems improbable to me. Lets check the smooths. Or only go to 1000m.
, legend.background = element_blank(), legend.box.background = element_rect()
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_b <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_c <- function(x) {
xout <- gsub("1e", "10^{", format(x),fixed=TRUE)
xout <- gsub("{-0", "{-", xout,fixed=TRUE)
xout <- gsub("{+", "{", xout,fixed=TRUE)
xout <- gsub("{0", "{", xout,fixed=TRUE)
xout <- paste(xout,"}",sep="")
return(parse(text=xout))
}
scale_x_log10nice <- function(name=NULL,omag=seq(-10,20),...) {
breaks10 <- 10^omag
scale_x_log10(breaks=breaks10,labels=scientific_10_c(breaks10),...)
}
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10nice()+
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
#labs(x = "moo", y = "Depth (m)") +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm")) +
geom_segment(aes(y = 160, yend = 160, x = 20, xend = 500), color = "darkgreen", stroke = 0.5)+
geom_segment(aes(y = 850, yend = 850, x = 20, xend = 500), color = "darkblue", stroke = 0.5)#+ geom_hline(yintercept = 850, color = "darkblue")
Ignoring unknown parameters: strokeIgnoring unknown parameters: stroke
pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)
Removed 14 rows containing missing values (geom_point).
pltFlx
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse()+
#scale_x_log10() +
scale_x_log10(breaks = c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20)), limits = c(20, 180)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")+
geom_hline(yintercept = 160, color = "darkgreen")
pltFlxZoom
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")+
geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3
#plotly::ggplotly(plt1pos)
# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
#
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
#
# pgTop <- plot_grid(pltFlxNoLegend, pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""), label_size = 14)
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
#
# pgBoth
#
# ggsave("figures/FluxDeepDive.png")
# ggsave("figures/FluxDeepDive.svg")
Within panel drawing
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","B"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgTop
pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)
I don’t know whats going on below here
pgBottom <- pltDelta3 + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = NA, inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)
Removed 33 rows containing missing values (geom_point).
pgBoth
ggsave("figures/FluxDeepDive.png")
Saving 5 x 9 in image
ggsave("figures/FluxDeepDive.svg")
Saving 5 x 9 in image
# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
#
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
#
# pgTop <- plot_grid(pltFlxNoLegend + ylim(c(1000, 0)), pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3 + ylim(c(1000, 0)), pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""))
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
#
# pgBoth
#
# #ggsave("figures/FluxShallowDive.png")
# #ggsave("figures/FluxShallowDive.svg")
Test for day to day and hourly variability in delta flux
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05697 -3.227 0.00267 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.716 1.919 4.652 0.03562 *
s(Day) 1.931 1.995 6.788 0.00412 **
s(Hour) 1.555 2.000 3.247 0.02318 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.328 Deviance explained = 41.3%
GCV = 0.15991 Scale est. = 0.1363 n = 42
summary(DFG2)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06186 -2.972 0.00514 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.641 1.871 3.599 0.0798 .
s(Day) 1.827 1.970 3.034 0.0463 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.207 Deviance explained = 27.4%
GCV = 0.17984 Scale est. = 0.1607 n = 42
summary(DFG3)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06635 -2.771 0.00849 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.592 1.834 3.213 0.107
R-sq.(adj) = 0.0878 Deviance explained = 12.3%
GCV = 0.19707 Scale est. = 0.18491 n = 42
summary(DFG1)$r.sq - summary(DFG2)$r.sq
[1] 0.1204131
summary(DFG2)$r.sq - summary(DFG3)$r.sq
[1] 0.1193917
summary(DFG3)$r.sq
[1] 0.08782536
png(filename = “./figures/CombinedP2Info.png”, width = 10, height = 8, units = “in”, res = 200) StationInfoPlot() dev.off()
#plot.new()
FluxGamPlot <- function(){
par(mfrow = c(2,2))
plot(DFG1)
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}
FluxGamPlot()
png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)
FluxGamPlot()
dev.off()
png
2
Check of actual data for hour
ggplot(data = bds %>% filter(depth >= 175, depth <= 500), aes(y = DF/DZ, x = hour(time), col = depth, group = depth)) + geom_point() + geom_line()
#Osps
(u mol C / m^3 / day)
disagFig <- bds %>% filter(project == "ETNP") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
#labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day)), y = "Depth (m)") +
labs(x = paste("OSMS", expression((μmol/m^3/day)))) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
disagFig
#ggsave("..figures/FluxSizeShift.svg"
ggsave("figures/FluxSizeShift.png")
Saving 6 x 4 in image
ggsave("figures/FluxSizeShift.svg")
Saving 6 x 4 in image
For NSF Proposal plot_grid( PlotParticlesmany, PlotPSDmany, rel_widths = c(2, 3) )
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","D"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgCF <- plot_grid(PlotParticlesmany + ylim(c(1000, 0)),
PlotPSDmany + ylim(c(1000, 0)) +theme(legend.position = "none"),
pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
disagFig + theme(plot.margin = unit(c(0, 0, .2, 0), units = "cm")),
ncol = 2, rel_heights = c(4, 4), labels = c("A", "B", "C", "E"), label_size = 16)
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Removed 308 rows containing missing values (geom_point).Removed 308 rows containing missing values (geom_point).Removed 37 rows containing missing values (geom_point).
pgCF
ggsave("./figures/UVP_for_zoop_proposal.png", pgCF, width = 10, height = 8, units = "in")
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(OZG1)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.023796 0.002244 10.6 6.42e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.9032 1.991 8.166 0.000545 ***
s(Day) 1.8673 1.982 4.932 0.015920 *
s(Hour) 0.1508 2.000 0.080 0.345843
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.304 Deviance explained = 35.1%
GCV = 0.00032912 Scale est. = 0.00030213 n = 60
summary(OZG2)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.023796 0.002247 10.59 6.54e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.904 1.991 8.157 0.000549 ***
s(Day) 1.868 1.983 4.868 0.016638 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.302 Deviance explained = 34.7%
GCV = 0.00032916 Scale est. = 0.00030298 n = 60
summary(OZG3)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0238 0.0024 9.916 5.03e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.894 1.989 7.458 0.000924 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.204 Deviance explained = 23%
GCV = 0.00036308 Scale est. = 0.00034556 n = 60
summary(OZG1)$r.sq - summary(OZG2)$r.sq
[1] 0.001957736
summary(OZG2)$r.sq - summary(OZG3)$r.sq
[1] 0.09804936
summary(OZG3)$r.sq
[1] 0.2043621
OSMSGamPlot <- function(){
par(mfrow = c(1,2))
plot(OZG2)
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
}
OSMSGamPlot()
png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 6, units = "in", res = 200)
OSMSGamPlot()
dev.off()
png
2
plot(OZG2)
bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = ospsDZ)) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(500, 0)) + geom_vline(xintercept = 0)
trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")
Parsed with column specification:
cols(
Class = [31mcol_character()[39m,
Depth = [32mcol_double()[39m,
TrapID = [31mcol_character()[39m,
TrapType = [31mcol_character()[39m,
SampleType = [31mcol_character()[39m,
C_flux = [32mcol_double()[39m,
C_flux_umol = [32mcol_double()[39m
)
UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
Flux = [32mcol_double()[39m
)
UVPFluxOE <- read_csv("dataOut/ObservedVsExpectedFlux.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
tn_flux = [32mcol_double()[39m,
profile = [31mcol_character()[39m,
project = [31mcol_character()[39m,
time = [31mcol_character()[39m,
tot_flux2 = [32mcol_double()[39m
)
trapFlux3
UVPFluxComb
fluxMS_distilled_toPlot <- trapFlux3 %>%
mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))
Remove traps where mass spec didn’t work correctly 2-17 150 1-12 73m 1-12 148 2-14 100 |(TrapID == “2-17” & Depth == 150)
fluxMS_distilled_toPlot2 <- fluxMS_distilled_toPlot %>%
filter(!((TrapID == "1-12") | (TrapID == "2-14" & Depth == 100)|(TrapID == "2-17" & Depth == 150)))
fluxMS_distilled_toPlot2
UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, fill = SampleType, shape = TrapType),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot) +
geom_point(aes(x = Flux), size = 3, shape = 21, color = "white", fill = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP")) + # dummy point for the legend
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") + xlab("Flux µmolC/m^2/day") +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
scale_fill_viridis_d() +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
) +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
Traps where mass spec didn’t work.
UVPFluxPlot00 <- UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, shape = TrapType, ID = TrapID),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot2) +
geom_line(aes(x = Flux), size = 1, color = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP Estimate")) + # dummy point for the legend
geom_point(aes(x = tot_flux2), size = 3, shape = 21, color = "white", fill = "black", data = UVPFluxOE) +
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") +
#xlab(expression(Flux µmolC/m^2/day)) +
xlab(expression(paste("x axis ", ring(A)^2))) +
xlab(expression(paste("Flux (µmolC/", m^2, "/day)"))) +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
)
Ignoring unknown aesthetics: ID
# UVPFluxPlot <- UVPFluxPlot00 +
# geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
UVPFluxPlot00
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
horizontalGamPlot <- dataGamHorizontal %>% ggplot(aes(x = resp_fit, y = depth, col = log(lb), group = lb)) + scale_y_reverse() + geom_point() + scale_x_log10(limits = c(10^-8, NA)) + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + geom_errorbar(aes(xmin = resp_lower, xmax = resp_upper), width = 10, alpha = 0.5)+ theme_bw()
TPPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = TotalParticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "TotalParticles Observed (#)")
nnpPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = n_nparticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Binsize and Volume Normalized Particles (#/L/mm)")
FitPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = nnp_smooth, xmin = nnp_lower, xmax = nnp_upper, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Smoothed - Normalized Particles (#/L/mm)") + geom_errorbar(width = 10, alpha = 0.5)
npLegend <- get_legend(FitPlot + theme(legend.box.margin = margin(0, 0, 40, 200)) + labs(col = expression(log[e](Size (mm)))))
Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
plot_grid(
TPPlot + theme(legend.position = "none"),
nnpPlot + theme(legend.position = "none"),
npLegend ,
FitPlot + theme(legend.position = "none")
)
Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
ggsave("figures/AllParticleSizes.svg")
Saving 10 x 6.18 in image
ggsave("figures/AllParticleSizes.png")
Saving 10 x 6.18 in image
SameGam <- gam(TotalParticles ~s(log(lb), log(depth), by = factor(time)), offset = log(vol * binsize), family = nb(),
data = bes %>% filter(project == "ETNP"))
besE <- bes %>% filter(project == "ETNP")
lb_new <- exp(seq(from = log(0.1), to = log(2.1), by = 0.05))
ub_new <- lead(lb_new)
binsize_new <- ub_new - lb_new
lbbs <- tibble(lb = lb_new, ub = ub_new, binsize = binsize_new)
Expanded <- expand_grid(lb = exp(seq(from = log(0.1), to = log(2), by = 0.05)), depth = seq(from = 20, to = 2000, by = 20), time = as.factor(unique(besE$time))) %>% left_join(lbbs, by = "lb")
Pred <- exp(predict(SameGam, Expanded))
ToPlot <- bind_cols(Expanded, nnparticles = Pred) %>% mutate(time = as.character(time)) %>% mutate(nparticles = nnparticles * binsize)
ToPlot %>% filter(lb <= 2) %>% ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c() + scale_y_reverse() + scale_x_log10() + facet_wrap(~time) + geom_contour(color = "black")
meanBese <- ToPlot %>% filter(lb <= 2) %>% group_by(lb, depth) %>% summarize(nparticles = mean(nparticles), nnparticles = mean(nnparticles))
WBColorMap <- meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = "log10(number density \n (normalized))") + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
WBColorMap
Average of everything
#meanBese043 <- ToPlot %>% filter(lb <= 2, time == "2017-01-13 11:51:31")
meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c() + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen")
Just 043
mbGam <- meanBese %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
mbGam <- meanBese043 %>% group_by(depth) %>% nest() %>% mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>% mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
pWBPSD <- mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
pWBPSD
bds %>% filter(profile == “stn_043”, depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = “darkgreen”) + geom_hline(yintercept = 950, color = “darkblue”) + geom_ribbon(alpha = 0.2) + labs(x = “PSD slope”)
All of them
bds %>% filter(profile == "stn_043", depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_ribbon(alpha = 0.2) + labs(x = "PSD slope")
043 only
bds %>% filter(profile == "stn_043", depth <= 2000, depth > 175) %>% ggplot(aes(x = small_biovolume, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_point()
ubDf0 <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global)
ubDf <- ubDf0 %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth)
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
ubDf %>% ggplot(aes(x = nbiomass, y = depth , group = time, col = time)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1))
ubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
ubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 175, color = "darkgreen")
PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- PubDf %>% filter(depth <= 165, depth >= 155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 160, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 850, color = "darkblue") + labs( x = "Small particle mass (norm.)")
pWBS
LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 165, depth >=155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 160, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 850, color = "darkblue")
pWBL
For tom and danielle
WBColorMap
pWBPSD
pWBS
pWBL
WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))
WBFig5
WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined
ggsave("figures/WBModelValidation.png")
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxP16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10(limits = c(35, 150),breaks = seq(from = 20, to = 150, by = 20)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
geom_hline(yintercept = 200, color = "darkgreen") +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))
#
#
#
# pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
# pltFlxLegend <- get_legend(pltFlx)
#
pltFlxP16
# #plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3P16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-1, .1), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 200, color = "darkgreen")+
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3P16
#plotly::ggplotly(plt1pos)
osms_p16 <- bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), group = factor(time))) + geom_point(size = 3) + geom_path() + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = "Observed - Modeled Small Particle Flux \n µmol/m^3/day") +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen")
plotly::ggplotly(osms_p16)
#ggsave("..figures/FluxSizeShift.svg"
plot_grid(
pltFlxP16,
pltDelta3P16,
osms_p16
)
ggsave("figures/P16FluxRelate.svg")
ggsave("figures/P16FluxRelate.png")
source("ModelStuff.R")
scan_for_example <- bds %>% filter(project == "ETNP", depth < 500, depth > 200) %>% select(profile, depth, DFP, use_this_DFP, ospsDZ)
loc_station = "stn_036"
loc_depth = 225
loc_prev_depth = loc_depth - 50
loc_DFP <- bds %>% filter(profile == loc_station, depth %in% c(275)) %>% pull(DFP)
loc_use_DFP <- bds %>% filter(profile == loc_station, depth %in% c(275)) %>% pull(use_this_DFP)
for_single_disag <- bes %>% filter(profile == loc_station, depth %in% c(225, 275)) %>% select(depth, lb, nnp_smooth) %>%
mutate(depth = recode(depth, `225` = "Shallow", `275` = "Deep")) %>%
pivot_wider(names_from = depth, values_from = nnp_smooth)
with_disag <- for_single_disag %>%
mutate(Predicted_Deep = remin_smooth_shuffle(Shallow, loc_use_DFP))
#remin_smooth_shuffle(for_single_disag$Shallow,loc_use_DFP)
for_plot_disag <- with_disag %>% pivot_longer(cols = -lb) %>% filter(lb <= 0.5) %>%
mutate(name = factor(name, levels = c("Shallow", "Deep", "Predicted_Deep"))) %>%
mutate(name = recode_factor(name, Shallow = "Shallow (250m)", Deep = "Deep (275m)", Predicted_Deep = "Predicted Deep (275m)"))
for_plot_disag %>% ggplot(aes(x = lb, y = value, shape = name)) + geom_point() + scale_x_log10() + scale_y_log10() + scale_shape_manual(values = c(1, 6, 3)) + theme(legend.title = element_blank()) + labs(x = "Particle Size (mm)", y = "Normalized Particle Abundance (#/L/mm)")
ggsave("figures/DisagExample.png")
ggsave("figures/DisagExample.svg")
dataBinned <- read_csv("data/backscatter_table_go7.csv")
dataBinned_01 <- dataBinned %>%
mutate(timeMex = with_tz(time_bin, tzone = "US/Central") )
startDay <- dataBinned_01$timeMex %>% na.omit %>% min %>% floor_date(unit = "days")
endDay <- dataBinned_01$timeMex %>% na.omit %>% max %>% ceiling_date(unit = "days")
timeBreaks <- seq(from = startDay, to = endDay, by = "12 hours")
timeLabels <- format(timeBreaks)
plotLetters <- tribble(
~letter, ~depth_bin, ~timeMex,
"A", 350, as.POSIXct("2017-01-07 13:00:00"),
"B", 200, as.POSIXct("2017-01-09 23:00:00"),
"C", 150, as.POSIXct("2017-01-08 13:00:00"),
"D", 625, as.POSIXct("2017-01-12 07:00:00"),
"E", 750, as.POSIXct("2017-01-10 02:00:00")
)
library(shadowtext)
plot18k <- dataBinned_01 %>% filter(frequency == 18000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue") +
geom_shadowtext(data = plotLetters, aes(x = timeMex - 12 * 60^2, y = depth_bin, label = letter), inherit.aes = FALSE, size = 6, bg.color = "white", color = "black") +
geom_segment(data = plotLetters, inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")), color = "white", size = 1.5) +
geom_segment(data = plotLetters, inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")))
plot18k
#ggsave("figures/stationP2_EK60_18kOnly.png")
plotLetters <- tribble(
~letter, ~depth_bin, ~timeMex,
"A", 300, as.POSIXct("2017-01-07 13:00:00"),
"B", 200, as.POSIXct("2017-01-09 23:00:00"),
"C", 150, as.POSIXct("2017-01-08 13:00:00"),
"D", 625, as.POSIXct("2017-01-12 07:00:00"),
"E", 750, as.POSIXct("2017-01-10 02:00:00")
)
library(shadowtext)
plot38k <- dataBinned_01 %>% filter(frequency == 38000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -95), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "Day::Hour", y = "Depth (m)", fill = "Backscatter (dB)") + theme_bw() + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue") +
geom_shadowtext(data = plotLetters[1:2,], aes(x = timeMex - 12 * 60^2, y = depth_bin, label = letter), inherit.aes = FALSE, size = 6, bg.color = "white", color = "black") +
geom_segment(data = plotLetters[1:2,], inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")), color = "white", size = 1.5) +
geom_segment(data = plotLetters[1:2,], inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc"))) +
theme(axis.text.x = element_text(size = 12,angle = 90, vjust = 0.5),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 14),
legend.title = element_text(size = 14)
)
plot38k
ggsave("figures/stationP2_EK60_38kOnly.png")
plot18k + scale_y_reverse(limits = c(500, 200))
plot200k <- dataBinned_01 %>% filter(frequency == 200000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -155), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
plot200k
dataBinned_01 %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "Depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
facet_wrap(~frequency, ncol = 1)
#ggsave("figures/stationP2_EK60_go7.svg", width = 4, height = 10)
ggsave("figures/stationP2_EK60_go7.png", width = 6, height = 10)